import numpy as np
import scipy as sp
import scipy.stats
import scipy.special
import scipy.integrate
import statsmodels as sm
from statsmodels.distributions.empirical_distribution import ECDF
import plotly
import plotly.graph_objs as go
from numba import njit
np.set_printoptions(
precision=5,
linewidth=120,
)
plotly.offline.init_notebook_mode()
Гипергеометрическое распределение - дискретное распределение, описывающее вероятность события, при котором ровно $k$ из $n$ случайно выбранных элементов окажутся помеченными, при этом выборка осуществляется из множества мощности $N$, в котором присутствует $m$ помеченных элементов. Считается, что каждый из элементов может быть выбран с одинаковой вероятностью $\frac{1}{N}$. Запишем это формально: $$\begin{gathered} N \in \mathbb{N},\ m \in \overline{0, N},\ n \in \overline{0, N},\\ k \in \overline{\max(0, m + n - N), \min(m, n)} \end{gathered}$$
Тогда $HG(N, m, n)$ описывает вероятность события, при котором ровно $k$ из $n$ элементов выборки окажутся помеченными: $$\left\{\ \xi \sim HG(N, m, n)\ \right\} \iff \left\{\mathbb{P}\left(\, \xi=k \,\right) = \frac{\binom{m}{k}\binom{N-m}{n-k}}{\binom{N}{n}}\right\}$$
По определению, математическое ожидание случайной величины – это ее $1^\text{й}$ начальный момент. Для начала, найдем $k^\text{й}$ начальный момент для $\xi$ (это понадобится для дальнейших выводов): $$\mathbb{E}\left[\, \xi^r \,\right] = \sum_{k=0}^{n} k^r \cdot \mathbb{P}\left(\, \xi=k \,\right) = \sum_{k=0}^{n} k^r\frac{\binom{m}{k}\binom{N-m}{n-k}}{\binom{N}{n}}$$ Можем считать, что сумма берется при $k=\overline{1,n}$, так как слагаемое при $k=0$ будет равно $0$. Заметим, что $$\begin{aligned} k\binom{m}{k} &= k \frac{m!}{k!(m-k)!} =\\ &= k \frac{m \cdot (m-1)!}{k \cdot (k-1)! \cdot (m-k)!} =\\ &= m \frac{(m-1)!}{(k-1)! \cdot \left(m-1 - (k-1)\right)!} =\\ &= m \binom{m-1}{k-1} \end{aligned}$$ и, как следствие, $$\binom{N}{n} = \frac{1}{n} \cdot n \cdot \binom{N}{n} = \frac{1}{n} N \binom{N-1}{n-1}$$ Подставим: $$\mathbb{E}\left[\, \xi^r \,\right] = \frac{n \cdot m}{N} \sum_{k=1}^{r-1} \frac{\binom{m-1}{k-1}\binom{N-m}{n-k}}{\binom{N-1}{n-1}}$$ Положим $j := k-1$ и изменим индекс суммирования на $j = \overline{0, n-1}$. Заметим, что $n - k = n - (j+1) = (n-1) - j$ и $N - m = (N-1) - (m-1)$: $$\mathbb{E}\left[\, \xi^r \,\right] = \frac{n \cdot m}{N} \color{dodgerblue}{\sum_{j=0}^{n-1} (j+1)^{r-1} \frac{\binom{m-1}{j}\binom{(N-1) - (m-1)}{(n-1) - j}}{\binom{N-1}{n-1}}}$$ Заметим, что выделенная часть выражения может быть записана, как $\mathbb{E}\left[\, (\theta+1)^{r-1} \,\right]$, где $\theta \sim HG(N-1, m-1, n-1)$. Следовательно, $$\mathbb{E}\left[\, \xi^r \,\right] = \frac{n \cdot m}{N} \mathbb{E}\left[\, (\theta+1)^{r-1} \,\right]$$ Таким образом, $$\boxed{ \mathbb{E}\left[\, \xi \,\right] = \frac{n \cdot m}{N} }$$
По определению дисперсии, $$\mathrm{Var}\left[\, \xi \,\right] = \mathbb{E}\left[\, \left(\ \xi - \mathbb{E}\left[\, \xi \,\right]\ \right)^2 \,\right] = \mathbb{E}\left[\, \xi^2 \,\right] - \left(\ \mathbb{E}\left[\, \xi \,\right]\ \right)^2$$ Выведем $2^\text{й}$ начальный момент: $$\mathbb{E}\left[\, \xi^2 \,\right] = \frac{n \cdot m}{N}\mathbb{E}\left[\, \theta+1 \,\right] = \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1}+1\right)$$ Подставим: $$\begin{aligned} \mathrm{Var}\left[\, \xi \,\right] &= \mathbb{E}\left[\, \xi^2 \,\right] - \left(\mathbb{E}\left[\, \xi \,\right]\right)^2 =\\ &= \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1}+1\right) - \left(\frac{n \cdot m}{N}\right)^2=\\ &= \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1} + 1 - \frac{n \cdot m}{N}\right) \end{aligned}$$ Таким образом, $$\boxed{ \mathrm{Var}\left[\, \xi \,\right] = \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1} + 1 - \frac{n \cdot m}{N}\right) }$$
По определению, производящая функция вероятностей $G(z, \xi)$ – это математическое ожидание новой случайной величины $z^\xi$. То есть: $$G_\xi(z) = \mathbb{E}\left[\, z^\xi \,\right]$$ Для $\xi \sim HG(N, m, n)$ производящая функция выглядит так: $$G_\xi(z) = \binom{N-m}{n}\left({}_{2}F_{1}(-m, -n; N-m-n+1; z) - 1\right)$$ Здесь ${}_{2}F_{1}$ - это гипергеометрическая функция, определенная следующим образом: $${}_{2}F_{1}(a,b;c;z) = \sum_{n=0}^{\infty} \frac{a^{(n)} b^{(n)}}{c^{(n)}} \frac{z^n}{n!}$$ , а $x^{(n)}$ - возрастающий факториал, определенный как: $$x^{(n)} = \prod_{k=0}^{n-1} (x + k)$$
По определению, характеристическая функция случайно величины $\xi$ задается следующим образом: $$\varphi_\xi (t) = \mathbb{E}\left[\, e^{it\xi} \,\right]$$ Для $\xi \sim HG(N, m, n)$ характеристическая функция выглядит так: $$M_\xi(t) = \frac{\binom{N-D}{n}}{\binom{N}{n}}\ {}_{2}F_{1}\left(-n, -D; N-D-n+1; e^{it}\right)$$ Здесь ${}_{2}F_{1}$ - это гипергеометрическая функция.
Гистограмма – это графическое представление функции, приближающей плотность вероятности распределения на основе выборки из него.
Чтобы построить гистограмму, сначала нужно разбить множество значений выборки на несколько отрезков. Чаще всего, берут отрезки одинаковой длины, чтобы облегчить восприятие получившегося результата, однако это необязательно. Далее подсчитывается количество вхождений элементов выборки в каждый из отрезков и рисуются прямоугольники, по площади пропорциональные количеству попавших элементов выборки в соответствующий отрезок.
Вообще говоря, гистограмму можно использовать не только для приближения плотности на основе выборки, но и для визуализации самой плотности распределения, зная его плотность.
Мы будем строить гистограмму вероятностей, писать будем на языке Python3 и использовать следующие библиотеки:
Итак, для начала, определим класс гипергеометрического распределения
HG, который будет содержать в себе информацию о параметрах $N$, $m$ и
$n$ и предоставлять метод p(k), возвращающий вероятность принятия
случайной величиной значения k при данных параметрах:
class HG(object):
def __init__(self, N: int, m: int, n: int):
self.N = N
self.m = m
self.n = n
def p(self, k: int) -> float:
return sp.special.comb(self.m, k) \
* sp.special.comb(self.N-self.m, self.n-k) \
/ sp.special.comb(self.N, self.n)
@property
def expected(self):
return self.n * self.m / self.N
@property
def variance(self):
return self.expected * ((self.n - 1) * (self.m - 1) / (self.N - 1) + 1 - self.expected)
@property
def domain(self):
return (max(0, self.m + self.n - self.N), min(self.m, self.n))
def __str__(self) -> str:
return f'HG({self.N}, {self.m}, {self.n})'
Далее создадим объект случайной величины $\xi \sim HG(30, 15, 20)$:
xi = HG(80, 15, 40)
Следующим шагом, определим интервал $\overline{0, n}$, на котором мы будем рисовать нашу гистограмму:
hg_data_x = np.arange(*xi.domain)
И, наконец, построим гистограмму и выведем ее:
theoretical_hg_PDF = np.apply_along_axis(xi.p, 0, hg_data_x,)
theoretical_hg_PDF_trace=go.Scatter(
name='PDF',
legendgroup='PDF',
x=hg_data_x,
y=theoretical_hg_PDF,
mode='markers',
marker_color='blue',
)
hg_hist_fig = go.Figure(
data=(theoretical_hg_PDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\xi \sim ' + str(xi) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$\mathbb{P}(\xi=k)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$k$',
),
),
),
)
plotly.offline.iplot(hg_hist_fig)
По определению, функция распределения $F_\xi(k) = \mathbb{P}\left(\, \xi < k \,\right)$. Для дискретной случайной величины событие $\{\xi < k\} = \bigcup\limits_{i=0}^{k-1}\{\xi=i\}$. Каждое из событий $\{\xi=i\}\; \forall i \in \overline{0, k-1}$ являются попарно несовместными. То есть $\forall i,j \in \overline{0, k-1}: i \neq j$ выполняется $\{\xi=i\}\cap\{\xi=j\}=\emptyset$. Из этого следует, что $$\mathbb{P}\left(\, \xi < k \,\right) = \sum_{i=0}^{k-1}\mathbb{P}\left(\, \xi = i \,\right)$$ Подставим и получим: $$F_\xi(k) = \sum_{i=0}^{k-1}\mathbb{P}\left(\, \xi = i \,\right) = \sum_{i=0}^{k-1}\frac{\binom{m}{i}\binom{N-m}{n-i}}{\binom{N}{n}}$$
Построим график этой функции, учитывая, что аргументом $k$ должно быть натуральное число, не превосходящее $n$:
theoretical_hg_CDF = np.cumsum(np.apply_along_axis(xi.p, 0, hg_data_x))
theoretical_hg_CDF_trace=go.Scatter(
name='CDF',
legendgroup='CDF',
x=hg_data_x,
y=theoretical_hg_CDF,
line=go.scatter.Line(
shape='hv',
color='blue',
),
)
hg_dist_fig = go.Figure(
data=(theoretical_hg_CDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\xi \sim ' + str(xi) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$\mathbb{P}(\xi<k)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$k$',
),
),
),
)
hg_dist_fig.show()
Нормальное распределение - непрерывное распределение, описывающее поведение величины отклонения измеряемого значения $x$ от истинного значения $\mu$ (которое является математическим ожиданием) и в рамках некоторого разброса $\sigma$ (среднеквадратичного отклонения). Запишем это формально: $$\left\{ \eta \sim N(\mu, \sigma^2) \right\} \iff \left\{\begin{gathered} F_\eta(x) = \mathbb{P}\left(\, \eta < x \,\right) = \int_{-\infty}^{x} f_\eta(x)dx,\\ \text{где } f_\eta(x) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} \text{– плотность вероятности} \end{gathered}\right\}$$
Найдем математическое ожидание $\eta \sim N(\mu, \sigma^2)$: $$\begin{aligned} \mathbb{E}\left[\, \eta \,\right] &= \int_{-\infty}^{+\infty} x \cdot f_\eta(x)dx =\\ &= \int_{-\infty}^{+\infty} xe^{-\frac{(x-\mu)^2}{2\sigma^2}}dx =\\ &= \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{+\infty} xe^{-\frac{(x-\mu)^2}{2\sigma^2}}dx \end{aligned}$$ Сделаем замену $t = \frac{x-\mu}{\sqrt{2}\sigma}$: $$\begin{aligned} \mathbb{E}\left[\, \eta \,\right] &= \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{+\infty}(\sigma\sqrt{2}t + \mu) e^{-t^2} d\left(\frac{x-\mu}{\sqrt{2}\sigma}\right) =\\ &= \frac{\sigma\sqrt{2}}{\sqrt{\pi}}\int_{-\infty}^{+\infty}te^{-t^2}dt + \frac{\mu}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-t^2}dt =\\ &= \frac{\sigma\sqrt{2}}{\sqrt{\pi}}\left(\int_{-\infty}^{0}te^{-t^2}dt - \int_{-\infty}^{0}te^{-t^2}dt\right) + \frac{\mu}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-t^2}dt =\\ &= \frac{\mu}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-t^2}dt \end{aligned}$$ Заметим, что получившееся выражение содержит интеграл, который может быть сведен к интегралу Эйлера-Пуассона: $$\int_{-\infty}^{+\infty}e^{-t^2}dt = 2\int_{0}^{+\infty}e^{-t^2}dt = \sqrt{\pi}$$ Таким образом, $$\boxed{ \mathbb{E}\left[\, \eta \,\right] = \mu }$$
Сделаем ту же замену переменной $t = \frac{x-\mu}{\sqrt{2}\sigma}$, тогда $x = t\sqrt{2}\sigma+\mu$ и: $$\begin{aligned} \mathrm{Var}\left[\, \eta \,\right] &= \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{+\infty}(\sqrt{2}\sigma)^2 t^2 e^{-t^2}d(t\sqrt{2}\sigma+\mu) =\\ &= \frac{2\sigma^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}t^2 e^{-t^2}dt \end{aligned}$$ Проинтегрируем по частям: $$\begin{aligned} \mathrm{Var}\left[\, \eta \,\right] &= \frac{\sigma^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}t 2t e^{-t^2} dt =\\ &= \frac{\sigma^2}{\sqrt{\pi}}\left(\left. -t e^{-t^2} \right|_{-\infty}^{+\infty} + \int_{-\infty}^{+\infty}e^{-t^2}dt\right) \end{aligned}$$ Здесь снова появляется интеграл Эйлера-Пуассона и, в итоге, получаем: $$\boxed{ \mathrm{Var}\left[\, \eta \,\right] = \sigma^2 }$$ То есть, $\sigma$ является среднеквадратичным отклонением.
Характеристическая функция для $\eta \sim N(\mu, \sigma^2)$ имеет вид: $$\varphi_\eta(t) = \exp\left(\mu it + \frac{\sigma^2 t^2}{2}\right)$$
Определим класс нормального распределения N, который будет содержать в
себе информацию о параметрах $\mu$ и $\sigma$ и предоставлять следующие
методы:
f(x) - возвращает значение плотности в точке xp(k) - возвращает
$\mathbb{P}\left(\, \eta < x \,\right) = \int_{-\infty}^x f_\eta(x)dx$class Norm(object):
def __init__(self, mu: float, sigma: float):
self.mu = mu
self.sigma = sigma
def f(self, x: float) -> float:
return np.exp(-((x-self.mu)**2)/(2*self.sigma**2))/(self.sigma*(2*np.pi)**.5)
def p(self, x: float) -> float:
return sp.integrate.quad(self.f, -np.inf, x)[0]
@property
def expected(self):
return self.mu
@property
def variance(self):
return self.sigma
def __str__(self):
return f'N({self.mu}, {self.sigma}^2)'
Далее создадим объект случайной величины $\xi \sim HG(30, 15, 20)$:
eta = Norm(0, 1)
Следующим шагом, руководствуясь правилом трех сигм, определим интервал $(-3\sigma, 3\sigma)$, в котором окажутся все значения случайной величины с вероятностью более $0.99$:
norm_data_x = np.linspace(-3*eta.sigma, 3*eta.sigma, 100)
И, наконец, построим плотность, используя метод N.f нашего класса:
theoretical_norm_PDF = np.vectorize(eta.f)(norm_data_x)
theoretical_norm_PDF_trace = go.Scatter(
name='PDF',
legendgroup='PDF',
x=norm_data_x,
y=theoretical_norm_PDF,
line=go.scatter.Line(
color='blue',
),
)
norm_dens_fig = go.Figure(
data=(theoretical_norm_PDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\eta \sim ' + str(eta) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$f_\eta(x)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$x$',
),
),
),
)
plotly.offline.iplot(norm_dens_fig)
По определению, функция распределения $F_\eta(x) = \mathbb{P}\left(\, \eta < x \,\right)$. Для непрерывной случайной она определяется как интеграл от функции плотности вероятности: $$\mathbb{P}\left(\, \eta < x \,\right) = \int_{-\infty}^x f_\eta (x) dx$$. Подставим и получим: $$F_\xi(k) = \int_{-\infty}^x \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} dx = \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^x e^{-\frac{(x-\mu)^2}{2\sigma^2}} dx$$
Построим график этой функции, используя метод Norm.p:
theoretical_norm_CDF = np.vectorize(eta.p)(norm_data_x)
theoretical_norm_CDF_trace = go.Scatter(
name='CDF',
legendgroup='CDF',
x=norm_data_x,
y=theoretical_norm_CDF,
line=go.scatter.Line(
color='blue',
),
)
norm_dist_fig = go.Figure(
data=(theoretical_norm_CDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\eta \sim ' + str(eta) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$\mathbb{P}(\eta<x)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$x$',
),
),
),
)
plotly.offline.iplot(norm_dist_fig)
Типичной интерпретацией гипергеометрического распределения является выборка без возвращения из множества элементов, некоторые из которых являются помеченными. Представим, что в нашем распоряжении имеется корзина, наполненная шарами двух цветов: черные и белые. Причём всего в корзине находится $N$ шаров, $m$ из которых – белые. Шары в корзине тщательно перемешиваются, чтобы каждый из них мог быть вытащен с одинаковой вероятностью $\frac{1}{N}$. Далее случайно вытаскиваются $n$ шаров без возвращения. Гипергеометрическое распределение описывает вероятность того, что среди вытащенных шаров ровно $k$ окажутся белыми.
Действительно, всего существует $\binom{N}{n}$ выборок размера $n$, $\binom{m}{k}$ способов выбрать $k$ помеченных объектов (белых шаров), и $\binom{N-m}{n-k}$ способов заполнить оставшиеся $n-k$ слотов непомеченными объектами (черными шарами). Таким образом, вероятность того, что среди $n$ вытащенных объектов окажется ровно $k$ помеченных, будет равна $\frac{\binom{m}{k}\binom{N-m}{n-k}}{\binom{N}{n}}$.
Случайная величина, имеющая гипергеометрическое распределение с параметром $n=1$ будет иметь распределение Бернулли с параметром $m/N$: $$\left\{\xi \sim HG(N,m,1)\right\} \implies \left\{\xi \sim B\left(\frac{m}{N}\right)\right\}$$ Действительно, при размере выборки равным единице, случайная величина, имеющая гипергеометрическое распределение, может принимать только два значения $k \in \{0, 1\}$. То есть, нам либо попадется помеченный элемент, либо нет. Тогда эти вероятности описывает распределение Бернулли с вероятностью успеха, равной отношению общего количества элементов к количеству помеченных.
Если зафиксировать размер выборки и количество помеченных элементов, а мощность множества, из которого ведется выборка, устремить к бесконечности, то гипергеометрическое распределение будет сходиться к биномиальному: $$HG(N, m, n) \underset{N\to\infty}{\longrightarrow} Bi\left(n, \frac{m}{N}\right)$$
Нормальное распределение описывает нормированную случайную величину, которая является суммой многих случайных слабо взаимосвязанных величин, каждая из которых вносит малый вклад относительно общей суммы. Это вытекает из центральной предельной теоремы.
Сумма двух независимых случайных величин, имеющих нормальное распределение, имеет распределение Коши: $$\begin{aligned} \xi\ &\sim\ N(\mu_1, {\sigma_1}^2)\\ \eta\ &\sim\ N(\mu_2, {\sigma_2}^2)\\ \xi+\eta\ &\sim\ C(\mu_1 + \mu_2, \sqrt{{\sigma_1}^2 + {\sigma_2}^2}) \end{aligned}$$
Сумма квадратов $k$ независимых стандартных нормальных случайных величин имеет распределение $\chi^2$ c $k$ степенями свободы: $$\begin{aligned} \forall i \in \overline{1, k}\quad \xi_i\ &\sim\ N(0, 1)\\ \sum_{i=1}^k \xi_i &\sim \chi^2(k) \end{aligned}$$
Натуральный логарифм логнормального распределения имеет нормальное распределение: $$\begin{aligned} \xi &\sim LogN(\mu,\sigma^2)\\ \ln\xi &\sim N(\mu,\sigma^2) \end{aligned}$$
Основные алгоритмы моделирования брались из книги Р. Н. Вадзинского "Справочник по вероятностным распределениям"
Объявим константы:
QS = (5, 10, 10**2, 10**3, 10**5)
N = 5
Добавим к нашему классу HG метод gen, который будет возвращать реализацию из $HG(N, m, n)$:
@njit
def gen_hg(N: int, m: int, n: int) -> int:
x = 0
for _ in range(n):
if np.random.rand() < m/N:
x += 1
m -= 1
N -= 1
return x
HG.gen = lambda self: gen_hg(self.N, self.m, self.n)
Поясню происходящее: x – это количество помеченных элементов из тех, что мы выбрали.
Изначально x = 0, потому что мы еще не выбрали ни одного элемента.
Далее начинаем выбирать n элементов. Каждый раз, вероятность того, что мы выберем помеченный элемент равна m/N. Мы используем функцию numpy.random.rand, которая возвращает действительное число в интервале $[0, 1)$, и проверяем, меньше ли значение, чем m/N. Если это значение меньше, это означает, что мы выбрали один из помеченных элементов. В этом случае, мы увеличиваем x на единицу и уменьшаем количество оставшихся помеченных элементов на 1, потому что выборка ведется без возвращения.
Также, вне зависимости от того, был ли выбранный объект помеченным, мы уменьшаем количество объектов, из которых ведется выборка, на единицу.
Этот процесс повторяется ровно n раз, так как всего должно быть выбрано n элементов.
Для ускорения работы функции (она будет часто вызываться при построении выборок), мы используем декоратор из прекрасной библиотеки numba, который при первом вызове функции скомпилирует нашу функцию в машинный код, и, при последующих вызовах, будет вызываться уже скомпилированная функция, которая выполняется во много раз быстрее (в нашем случае прирост в скорости будет около 2000%).
Объем выборки будем обозначать $q$, чтобы не путать с параметром распреления $n$.
Напишем новый метод HG.sample_hg(q), который будет возвращать выборку объема $q$ из $HG(N, m, n)$:
def sample_hg(N: int, m: int, n: int, q: int) -> np.ndarray:
return np.fromiter((gen_hg(N, m, n) for _ in range(q)), np.int)
HG.sample = lambda self, q: sample_hg(self.N, self.m, self.n, q)
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ сгенерируем по 5 выборок $X_q^i$ объема $q$, где $i\in\overline{1, 5}$
samples_hg = {}
for q in QS:
samples_hg[q] = np.ndarray(shape=(N, q), dtype=np.int)
for i in range(N):
samples_hg[q][i] = xi.sample(q)
Теперь samples_norm – словарь, где ключом является объемом выборки, а значением по соответствующему ключу – массив из N = 5 выборок соответствующего объема из распределения $HG(80, 15, 20)$.
Выведем выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(samples_hg[q]):
print(f'#{i+1}: {s}')
Для каждой выбоки вычислим эмпирическую функцию распределения (ECDF) с помощью функции statsmodels.distributions.empirical_distribution.ECDF:
HG_ECDFs = {}
for q in samples_hg:
HG_ECDFs[q] = np.ndarray(shape=(samples_hg[q].shape[0], hg_data_x.shape[0]), dtype=np.float)
for i, sample in enumerate(samples_hg[q]):
HG_ECDFs[q][i] = ECDF(sample)(hg_data_x)
Теперь для каждой выборки построим эмпирическую функцию распределения, а также сравним получившуюся функцию с теоретической функцией распределения:
HG_ECDFs_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes='all',
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$k$',
y_title=r'$\mathbb{P}\left(\,\xi\leq k\,\right)$')
for i, q in enumerate(HG_ECDFs):
for j, sample in enumerate(HG_ECDFs[q]):
theoretical_hg_CDF_trace.showlegend = (i + j == 0)
HG_ECDFs_fig.add_trace(theoretical_hg_CDF_trace, row=i+1, col=j+1)
HG_ECDFs_fig.add_trace(
go.Scatter(
name='ECDF',
x=hg_data_x,
y=HG_ECDFs[q][j],
line=go.scatter.Line(
shape='hv',
color='red',
),
legendgroup='CDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
HG_ECDFs_fig.update_layout(height=800, width=900)
HG_ECDFs_fig.show()
Все в порядке, ECDF стремится к CDF при увеличении объема выборки.
Для каждого $q\in\{5, 10, 100,1000, 10^5\}$ найдем верхнюю границу разности каждой пары эмпирических функций распределения $\Delta_q^{i,j}$: $$\Delta_q^{i,j} = \sup\limits_{y\,\in\,\rm{supp}\,\xi}\left(\,\hat{F}\left(y, \overrightarrow{x_i}\right) - \hat{F}\left(y, \overrightarrow{x_j}\right)\,\right), \quad i,j\,\in\overline{0, 5}$$
for q in HG_ECDFs:
print(f'====== q = {q} ======:')
for i, ECDF_i in enumerate(HG_ECDFs[q]):
for j, ECDF_j in enumerate(HG_ECDFs[q][i+1:], start=i+1):
print(f'i = {i}; j = {j}; max delta = {np.absolute(ECDF_i - ECDF_j).max():.4f}')
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим вариационный ряд выборки:
HG_VAR_ROWS = {}
for q in samples_hg:
HG_VAR_ROWS[q] = np.ndarray(shape=samples_hg[q].shape, dtype=np.int)
for i, sample in enumerate(samples_hg[q]):
HG_VAR_ROWS[q][i] = np.sort(sample)
Выведем получившиеся упорядоченные выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(HG_VAR_ROWS[q]):
print(f'#{i+1}: {s}')
Квантилью уровня $\alpha\in(0, 1)$ случайной величины $\xi$ называется такое число $x_\alpha\in\mathbb{R}$, что $$\begin{aligned} \mathbb{P}(\,\xi \leq x_\alpha\,) &\geq \alpha \\ \mathbb{P}(\,\xi \geq x_\alpha\,) &\geq 1-\alpha \\ \end{aligned}$$
Найдем квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ для гипергеометрического распределения $HG(80, 15, 40)$:
hg_quantiles = {}
for a in (.1, .5, .7):
for i, y in enumerate(theoretical_hg_CDF):
if y >= a:
hg_quantiles[a] = hg_data_x[i]
print(f'a = {a}, quantile = {hg_quantiles[a]}')
break
Выборочная квантиль уровня $\alpha\in(0, 1)$ выборки $\overrightarrow{X}$ - элемент вариационного ряда этой выборки стоящий на позиции $\left[\alpha\left|\overrightarrow{X}\right| + 1\right]$.
Для каждой выборки найдем квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ и сравним их с соответствующими квантилями данного распределения:
for q in QS:
print(f'======= q = {q} =======:')
for a in (.1, .5, .7):
print(f'------ a = {a} ------')
for i, sample in enumerate(HG_VAR_ROWS[q]):
quantile = sample[int(a*sample.shape[0])]
print(f'#{i}; quantile = {quantile:2}, real delta = {quantile - hg_quantiles[a]:+}')
Как и должно было быть, разность между выборочными квантилями и теоретическими стремится к нулю.
Вообще говоря, для вычисления квантилей выборок можно было воспользоваться функцией
np.quantile, но это было бы слишком просто...
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим гистограмму, а также сравним полученные графики с функцией вероятности $\mathbb{P}\left(\,\xi=k\,\right)$. Полигон частот строить не будем, потому что по сути, это просто альтернативный способ представления гистограммы, который будет нам только мешать анализировать графики.
HG_hists_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes=True,
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$k$',
y_title=r'$\mathbb{P}\left(\,\xi = k\,\right)$')
for i, q in enumerate(samples_hg):
for j, sample in enumerate(samples_hg[q]):
theoretical_hg_PDF_trace.showlegend = (i + j == 0)
HG_hists_fig.add_trace(theoretical_hg_PDF_trace, row=i+1, col=j+1)
HG_hists_fig.add_trace(
go.Histogram(
name='Histogram',
x=sample,
histnorm='probability',
marker_color='red',
legendgroup='PDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
HG_hists_fig.update_layout(height=800, width=900)
HG_hists_fig.show()
Как мы видим, закон больших чисел подтверждается графиками: с ростом объема выборки, значения гистограммы выборок все сильнее приближаются к истинному функции распределения.
С помощью формулы из справочника по вероятностным распределениям, мы можем получить сразу две реализации $x_i$ и $x_{i+1}$ из стандартного нормального распределения $N(0, 1)$, имея две реализации $r_i$ и $r_{i+1}$ из стандартного равномерного распределения на отрезке $[0, 1]$: $$\begin{aligned} x_i &= \sqrt{-2 \ln r_i} \sin(2\pi r_{i+1}) \\ x_{i+1} &= \sqrt{-2 \ln r_i} \cos(2\pi r_{i+1}) \end{aligned}$$
Чтобы получить реализацию $y_i$ из $N(\mu, \sigma)$: $$y_i = \mu + \sigma x_i$$
Для простоты, мы будем использовать только одну $y_i$ из двух возможных за раз:
@njit
def gen_norm(mu: float, sigma: float) -> float:
return mu + sigma * \
np.sqrt(-2 * np.log(np.random.rand())) * \
np.cos(2 * np.pi * np.random.rand())
Norm.gen = lambda self: gen_norm(self.mu, self.sigma)
Как говорилось ранее, объем выборки мы будем обозначать $q$.
Напишем метод HG.sample_norm(q), которsq будет возвращать выборку объема $q$ из $N(\mu, \sigma)$:
def sample_norm(mu: float, sigma: float, q: int) -> np.ndarray:
return np.fromiter((gen_norm(mu, sigma) for _ in range(q)), np.float)
Norm.sample = lambda self, q: sample_norm(self.mu, self.sigma, q)
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ сгенерируем по 5 выборок $X_\eta^{q,i}$ объема $q$, где $i\in\overline{1, 5}$
samples_norm = {}
for q in QS:
samples_norm[q] = np.ndarray(shape=(N, q), dtype=np.float64)
for i in range(N):
samples_norm[q][i] = eta.sample(q)
Теперь samples_norm – словарь, где ключом является объемом выборки, а значением по соответствующему ключу – массив из N = 5 выборок соответствующего объема.
Выведем выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(samples_norm[q]):
print(f'#{i+1}: {s}')
Для каждой выбоки вычислим эмпирическую функцию распределения (ECDF) с помощью функции statsmodels.distributions.empirical_distribution.ECDF:
Norm_ECDFs = {}
for q in samples_norm:
Norm_ECDFs[q] = np.ndarray(shape=(samples_norm[q].shape[0], norm_data_x.shape[0]), dtype=np.float)
for i, sample in enumerate(samples_norm[q]):
Norm_ECDFs[q][i] = ECDF(sample)(norm_data_x)
Теперь для каждой выборки построим эмпирическую функцию распределения, а также сравним получившуюся функцию с теоретической функцией распределения:
Norm_ECDFs_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes='all',
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$x$',
y_title=r'$\mathbb{P}\left(\,\eta\leq x\,\right)$')
for i, q in enumerate(Norm_ECDFs):
for j, sample in enumerate(Norm_ECDFs[q]):
theoretical_norm_CDF_trace.showlegend = (i + j == 0)
Norm_ECDFs_fig.add_trace(theoretical_norm_CDF_trace, row=i+1, col=j+1)
Norm_ECDFs_fig.add_trace(
go.Scatter(
name='ECDF',
x=norm_data_x,
y=Norm_ECDFs[q][j],
line=go.scatter.Line(
color='red',
),
legendgroup='CDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
Norm_ECDFs_fig.update_layout(height=800, width=900)
Norm_ECDFs_fig.show()
Тут опять все хорошо, значения ECDF стремятся к CDF при стремлении объема выборки к бесконечности.
Для каждого $q\in\{5, 10, 100,1000, 10^5\}$ найдем верхнюю границу разности каждой пары эмпирических функций распределения $\Delta_q^{i,j}$: $$\Delta_q^{i,j} = \sup\limits_{y\,\in\,\rm{supp}\,\eta}\left(\,\hat{F}\left(y, \overrightarrow{x_i}\right) - \hat{F}\left(y, \overrightarrow{x_j}\right)\,\right), \quad i,j\,\in\overline{0, 5}$$
for q in Norm_ECDFs:
print(f'====== q = {q} ======:')
for i, ECDF_i in enumerate(Norm_ECDFs[q]):
for j, ECDF_j in enumerate(Norm_ECDFs[q][i+1:], start=i+1):
print(f'i = {i}; j = {j}; max delta = {np.absolute(ECDF_i - ECDF_j).max():.4f}')
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим вариационный ряд выборки:
Norm_VAR_ROWS = {}
for q in samples_norm:
Norm_VAR_ROWS[q] = np.ndarray(shape=samples_norm[q].shape, dtype=np.float64)
for i, sample in enumerate(samples_norm[q]):
Norm_VAR_ROWS[q][i] = np.sort(sample)
Выведем получившиеся упорядоченные выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(Norm_VAR_ROWS[q]):
print(f'#{i+1}: {s}')
Найдем квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ для нормального распределения $N(0, 1^2)$:
norm_quantiles = {}
for a in (.1, .5, .7):
for i, y in enumerate(theoretical_norm_CDF):
if y >= a:
norm_quantiles[a] = norm_data_x[i]
print(f'a = {a}, quantile = {norm_quantiles[a]:+.4f}')
break
Теперь для каждой выборки найдем выборочные квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ и сравним их с соответствующими квантилями данного распределения:
for q in QS:
print(f'======= q = {q} =======:')
for a in (.1, .5, .7):
print(f'------ a = {a} ------')
for i, sample in enumerate(Norm_VAR_ROWS[q]):
quantile = sample[int(a*sample.shape[0])]
print(f'#{i}; quantile = {quantile:+2.3f}, real delta = {quantile - norm_quantiles[a]:+2.3f}')
Как и в случае с гипергеометрическим распределением, разность между выборочными квантилями и теоретическими стремится к нулю, что есть правильно.
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим гистограмму, а также сравним полученные графики с плотностью распределения $f_\eta(x)$. Как и в соответствующем пункте про гипергеометрическое распределение, полигон частот строить не будем.
Norm_hists_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes=True,
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$x$',
y_title=r'$f_\eta(x)$')
for i, q in enumerate(samples_norm):
for j, sample in enumerate(samples_norm[q]):
theoretical_norm_PDF_trace.showlegend = (i + j == 0)
Norm_hists_fig.add_trace(theoretical_norm_PDF_trace, row=i+1, col=j+1)
Norm_hists_fig.add_trace(
go.Histogram(
name='Histogram',
x=sample,
histnorm='probability density',
marker_color='red',
legendgroup='PDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
Norm_hists_fig.update_layout(height=800, width=900)
Norm_hists_fig.show()